function [im,inVector] = dct( I , mask)
%
step = 8;

if nargin <2
    mask = [1   1   1   1   1   0   0   0
            1   1   1   1   0   0   0   0
            1   1   1   0   0   0   0   0
            1   1   0   0   0   0   0   0
            1   0   0   0   0   0   0   0
            0   0   0   0   0   0   0   0
            0   0   0   0   0   0   0   0
            0   0   0   0   0   0   0   0];
end;

idx = find(mask ==1);
numOneBox = size(idx,1);


image = im2uint8(I);
d = size(image);
% zarovnat radky na nasobek step
if (ceil(d(1)/step) ~= d(1)/step)
    i = ceil(d(1)/step)* step - d(1);
    image = [image; zeros(i,d(2))];
end;
d = size(image);
% zarovnat sloupce na nasobek 8
if (ceil(d(2)/step) ~= d(2)/step)
    i = ceil(d(2)/step)* step - d(2);
    image = ([image zeros(d(1),i)]);
end;
d = size(image);
image = double(image) - 128;

sqrt1over8 = sqrt(1/8);
sqrt2over8 = sqrt(2/8);
piover8 = pi / 8;

Q = [16 11 10 16  24  40  51  61;
     12 12 14 19  26  58  60  55;
     14 13 16 24  40  57  69  56;
     14 17 22 29  51  87  80  62;
     18 22 37 56  68 109 103  77;
     24 35 55 64  81 104 113  92;
     49 64 78 87 103 121 120 101;
     72 92 95 98 112 100 103  99];

 x= repmat((0:step-1)',1,step);
 y = repmat(0:step-1,step,1);
 im = zeros(size(image));
 D = dctmtx(step);
 D_transposed = D';
 
inVector = zeros(d(1)/8 * d(2)/8 * numOneBox,1);
box = 1;
for j = 1:step:d(2)-1
    for i= 1:step:d(1)-1
        %j:j+step-1
        %i:i+step-1
        A = image(i:i+step-1,j:j+step-1);
        G = D * A * D_transposed;
        %G = zeros(size(A));
        %for u=0:1:7
        %    alphau = sqrt2over8;
        %    if u == 0
        %        alphau = sqrt1over8;
        %    end;
        %    for v=0:1:7
        %        alphav = sqrt2over8;
        %        if v == 0
        %            alphav = sqrt1over8;
        %        end;
        %        G(u+1,v+1) = alphau * alphav * sum(sum(A.*cos(piover8*(x + 1/2)*u) .* cos(piover8*(y + 1/2)*v) ));
        %    end;
        %end;
        DC = G ./ Q;
        im(i:i+step-1,j:j+step-1) = DC;
        inVector(box:box + numOneBox - 1,1) = DC(idx);
        box = box + numOneBox;
    end;
end;
    % viz anglicka wikipedie - jpeg
    % im = round(im);